sink("log/dbw_sim_script_log.txt", append = FALSE, type = "output")


## You can install the development version of dbw from GitHub:
devtools::install_github("hirotokatsumata/dbw")

## For simulation code
library(dbw)
library(furrr)
library(Hmisc)
library(tidyverse)


source("functions/dbw_simulation_functions.R")

## Set the number of cores
## Setting it as larger than 4 does not improve the speed
plan(strategy = multisession(workers = 4))

## Number of simulation iterations
B <- 2000

## Tables I.1--6 (Online Supplementary Materials I)
## Note: Outcome model is misspecified
result <- sim_main(B = B, save = TRUE, seed = 12345L, lambda200 = 0.03, lambda1000 = 0.007)

reslinear1 <- result$result$res200$lin1 %>%
              bind_rows(result$result$res1000$lin1)
reslinear2 <- result$result$res200$lin2 %>%
              bind_rows(result$result$res1000$lin2)
resquad1 <- result$result$res200$quad1 %>%
            bind_rows(result$result$res1000$quad1)
resquad2 <- result$result$res200$quad2 %>%
            bind_rows(result$result$res1000$quad2)
resexp1 <- result$result$res200$exp1 %>%
           bind_rows(result$result$res1000$exp1)
resexp2 <- result$result$res200$exp2 %>%
           bind_rows(result$result$res1000$exp2)

## Check convergence
result$notconverged


## Tables 2--4
sim_tbl2(result1 = reslinear1, result2 = reslinear2, 
         type = "Linear")
sim_tbl2(result1 = resquad1, result2 = resquad2, 
         type = "Quadratic")
sim_tbl2(result1 = resexp1, result2 = resexp2, 
         type = "Exponential")

result <- rbind(reslinear1, reslinear2, resquad1, resquad2, resexp1, resexp2) %>%
          filter(N == 1000 & 
                 Estimator %in% c("nDBW DR", "CBPS DR/BRDR", "CBPS DR", "Calibrated weighting DR", "Entropy balancing DR") & 
                 psmodel != "correct")
result$Estimator[result$Estimator == "nDBW DR"] <- "nDBW"
result$Estimator[result$Estimator == "CBPS DR/BRDR"] <- "CBPS"
result$Estimator[result$Estimator == "CBPS DR"] <- "CBPS"
result$Estimator[result$Estimator == "Calibrated weighting DR"] <- "CAL"
result$Estimator[result$Estimator == "Entropy balancing DR"] <- "EB"

## Figure 4 (left panel)
## Absolute bias
maxx <- max(abs(subset(result, Estimator == "nDBW", select = c("Bias", "RMSE")))) * 1.05
maxy <- max(abs(subset(result, Estimator %in% c("CAL", "CBPS", "EB"), select = c("Bias", "RMSE")))) * 1.05
result %>%
select(-c(RMSE)) %>%
pivot_wider(names_from = Estimator, values_from = Bias) %>%
pivot_longer(cols = c(CAL, CBPS, EB), names_to = "Estimator", values_to = "Bias") %>%
ggplot(., aes(x = abs(nDBW), y = abs(Bias), shape = Estimator, size = Estimator)) +
  geom_point(stroke = 0.75, alpha = 0.65) +
  geom_abline(intercept = 0, slope = 1) +
  xlab("nDBW DR estimator") +
  ylab("Calibrated weighting DR, CBPS DR, or entropy balancing DR estimator") +
  ggtitle("Absolute Bias") +
  scale_shape(name = "shape_size", solid = FALSE, 
              labels = c(CAL = "Calibrated weighting", CBPS = "CBPS", EB = "Entropy balancing")) +
  scale_size_manual(name = "shape_size", values = c(1.8, 1.5, 1.4), 
                    labels = c(CAL = "Calibrated weighting", CBPS = "CBPS", EB = "Entropy balancing")) +
  scale_x_continuous(expand = c(0, 0), limit = c(0, maxx)) +
  scale_y_continuous(expand = c(0, 0), limit = c(0, maxy)) +
  coord_fixed() +
  theme_classic() +
  theme(legend.position = c(0.28, 0.8),
        axis.title = element_text(size = 10.4),
        axis.text = element_text(size = 11),
        plot.title = element_text(hjust = 0.5, size = 12),
        legend.title = element_text(size = 0),
        legend.text = element_text(size = 10),
        legend.background = element_rect(fill = NA, colour = NA))
ggsave("figures/fig_4_left.pdf", width = 8.5, height = 12, units = "cm")
#ggsave("figures/sim_bias.pdf", width = 8.5, height = 12, units = "cm")

## Figure 4 (right panel)
## RMSE
maxx <- max(abs(subset(result, Estimator == "nDBW", select = c("Bias", "RMSE")))) * 1.05
maxy <- max(abs(subset(result, Estimator %in% c("CAL", "CBPS", "EB"), select = c("Bias", "RMSE")))) * 1.05
result %>%
select(-c(Bias)) %>%
pivot_wider(names_from = Estimator, values_from = RMSE) %>%
pivot_longer(cols = c(CAL, CBPS, EB), names_to = "Estimator", values_to = "RMSE") %>%
ggplot(., aes(x = abs(nDBW), y = abs(RMSE), shape = Estimator, size = Estimator)) +
  geom_point(stroke = 0.75, alpha = 0.65) +
  geom_abline(intercept = 0, slope = 1) +
  xlab("nDBW DR estimator") +
  ylab("Calibrated weighting DR, CBPS DR or entropy balancing DR estimator") +
  ggtitle("RMSE") +
  scale_shape(name = "shape_size", solid = FALSE, 
              labels = c(CAL = "Calibrated weighting", CBPS = "CBPS", EB = "Entropy balancing")) +
  scale_size_manual(name = "shape_size", values = c(1.8, 1.5, 1.4), 
                    labels = c(CAL = "Calibrated weighting", CBPS = "CBPS", EB = "Entropy balancing")) +
  scale_x_continuous(expand = c(0, 0), limit = c(0, maxx)) +
  scale_y_continuous(expand = c(0, 0), limit = c(0, maxy)) +
  coord_fixed() +
  theme_classic() +
  theme(legend.position = c(0.30, 0.8),
        axis.title = element_text(size = 10.4),
        axis.text = element_text(size = 11),
        plot.title = element_text(hjust = 0.5, size = 12),
        legend.title = element_text(size = 0),
        legend.text = element_text(size = 10),
        legend.background = element_rect(fill = NA, colour = NA))
ggsave("figures/fig_4_right.pdf", width = 8.5, height = 12, units = "cm")
#ggsave("figures/sim_rmse.pdf", width = 8.5, height = 12, units = "cm")


ratio <- result %>% mutate(id = paste0(type, posneg, psmodel)) %>%
select(c("id", "Estimator", "RMSE")) %>%
pivot_wider(id_cols = id,
            names_from = "Estimator",
            values_from = "RMSE") %>%
pivot_longer(cols = -c(id, nDBW)) %>%
mutate(ratio = value / nDBW) %>%
filter(name == "EB")

## Values mentioned in Section 5
ratio
mean(1 / ratio$ratio < 0.7)
mean(ratio$ratio < 0.7)

sink()
